Monte Carlo Simulations of Ion Transport Including Thermal Collisions: Collision 
Frequency Calculation and Benchmark Results 
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We have derived exact analytical formulas for piecewise calculation of collision frequency integrals 
for low energy ion transport in a background gas and in a mixture of background gases. These 
calculations are required for Monte Carlo simulations of ion transport and hybrid models of non 
equilibrium plasmas. The range of energies where it is necessary to apply this technique has been 
denned. Our calculations for the transport of ions in a model gas having a constant cross section are 
in excellent agreement with the benchmark data from accurate solutions to the Boltzmann equation. 
Benchmark calculations for ion motion in hot background argon are also given. 
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I. INTRODUCTION 

Non equilibrium plasmas are used in a broad range of 

f plications ranging from plasma etching [l|, Q, SI 0, d, 
0], plasma displays 
plications [H, 1 13i ] and 
Recently, microdischarges have become a popular tar- 
get for studies and development of new applications 
[18l [T9L I20I l21lj of non-equilibrium plasmas. Develop- 
ment of devices based on non-equilibrium plasmas is 
quite expensive and cannot be based on empirical meth- 
ods sopredictive accurate models had to be developed 

ssiii. 

Ion transport coefficients are used in both fluid and hy- 
brid models of plasmas @, SI SI SI, S3, SI, S3, S3, M, 
32l . |33| relevant for the applications mentioned above. In- 
directly transport coefficients are used to verify the com- 
pleteness and the absolute magnitude of the cross sec- 
tions in the sets [3, S3 that are to be used in Monte Carlo 
and Particle in Cell (PIC) models of plasmas [H, |35| . 
There is a critical shortage of data for most ions in most 
gases that are interesting for applications [361 ]. In addi- 
tion, there is a pressing need to have both accurate and 
efficient algorithms for simulation of ion motion. For ex- 
ample, even in hybrid models, that were developed to 
take into account the non-local transport of electrons, 
more and more frequently ions and even fast neutrals are 
treated by a Monte Carlo technique [H, S3, SI, S3 • 

In calculation of ion properties in a neutral back- 
ground gas with Monte Carlo simulation |40| or momen- 
tum transfer theory [4l[, it is very often assumed that 
the particles of the neutral gas are at rest. It is a reason- 
able approximation when the mean ion energy exceeds 
gas thermal energy by a large factor. This approxima- 



tion is, however, not valid for low electric fields, when the 
mean ion energy is comparable to the thermal energy of 
the background gas. In such a case collision frequency 
has to be calculated more accurately, based on relative 
velocity that takes into account the thermal motion of 
gas particles. 

Collision frequency for a test ion particle of velocity v 
in a neutral, background, gas is 



u(v) 



o-(|v-u|)|v-u|F(u)du, 



(1) 
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where n is the neutral gas density, F(u) is the gas velocity 
distribution function and a (|v — u|) is the cross section 
for ion-neutral scattering as a function of relative veloc- 
ity. Having in mind that collision frequency should be 
calculated for many different ion velocities v, one needs 
to be able to calculate ([1]) very efficiently. In addition, 
one should be aware that ions have thermal energies for a 
very wide range of reduced fields covering most practical 
discharges. 

Different approaches to calculation of the collision fre- 
quency can be found in the literature. In [42| calculation 
of the integral in (fTJ) is performed by numerical integra- 
tion using a Monte Carlo technique. While being poten- 
tially accurate, this technique is not a very efficient algo- 
rithm. Here, accuracy could be improved at the expense 
of the efficiency. Another approach is proposed in [43| 
where authors define the mean relative velocity, taking 
the cross section outside the integral in (TJ). Then they 
approximate the cross section by a cross section taken 
for that mean velocity. In this paper we have set out to 
find a more accurate yet still sufficiently efficient algo- 
rithm which may be used for Monte Carlo simulations of 
a quasi thermal ion (and electron) motion. In doing so 
we should remember that most cross sections are defined 
numerically at a limited number of points and that linear 
interpolation between those points is often used. Using 
the proposed algorithm we have also proposed bench- 
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mark calculations for ion motion in their parent gas at 
non-zero temperatures. 

After introduction, in section II we derive and analyse 
an explicit expression for collision frequency. Also we 
give an asymptotic formula for collision frequency which 
is valid for large ion velocities with respect to the back- 
ground gas thermal velocity. In section III we emphasise 
a proper way of sampling velocities of the background 
gas and propose rejection methods which depend on the 
cross section. Section IV deals with mixtures of gases. 
Numerical results, comparison with other data and pos- 
sible benchmark results are given in section V which is 
followed by conclusions. 



II. THEORY 

Explicit calculation of ([T|) is possible for specific cross 
sections and gas velocity distribution functions. We will 
assume that the background gas is in equilibrium at the 
temperature T, so it may be described by a Maxwellian 
velocity distribution function (MVDF). Introducing (the 

1/2 

most probable) gas thermal velocity w = (2kT/M) , 
where M is the mass of gas molecules and k is the Boltz- 
mann constant and after integration over solid angles, 
expression (fT]) can be transformed to 



v{v) 



p + OC 




(x — v) 2 






/ o a(x)x 2 


exp 
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exp 



(x + v) 2 



| Ax. (2) 



This expression is valid only for non-zero velocities of 
the ion. A similar expression can be derived for the ion 
velocity equal to zero, but this case is rarely encountered 
and will not be considered explicitly here. Alternatively, 
one can determine ^(0) as the limiting case v — > of {2]). 

It is possible to convert the integral of the type ((2]), to 
give the following expression for the collision frequency 
when ion velocities are much greater than the thermal 
gas velocity (see appendix for details): 



v(v) — na(v)v 1 + 



w 2 (v 2 a(v))' 



(3) 



We see that the first order correction is proportional to 
the square of the ratio w/v. Collision frequency for the 
cold gas is a special case w = of ((3|): 



v{v) = na(v)v. 



(4) 



This formula is very widely used in Monte Carlo simu- 
lations, for example in [4(| H, |45| even though it is not 
applicable at low ion energies. 

For further calculations of ([2]) we need a specific form 
of the cross section. Real cross sections are always tabu- 
lated as a function of incident ion energies (or velocities) 
or relative energies (velocities). In principle it could be 
argued that the cross sections are often tabulated as con- 
stant values in narrow energy bins so the solution for a 
constant cross section may be sufficient. However in the 
low energy limit the cross sections are often increasing 
rapidly and application of a constant cross section may 
require application of a large number of energy bins. At 
the same time usually the data for low energy limit of 
mobilities is available and thus one needs accurate and ef- 
ficient schemes for calculation in this energy range. Thus 
we assume that the cross section, given as a function of 
the relative ion-neutral velocity x, may be expressed as 



a(x) 



«?(!); 



-id) 
•4(f) 



0<x< Xl 



X\ < X < X2 



(5) 



«f (f) 2 +0* Xi-i<X< Xi 



Taking into account (JSJ) we can define integrals which 
originate in Q as: 



'1TWV 



2+m 



exp 



(x ± vf 



w 



dec, 
(6) 



with 77i == 0,1,2. These integrals may be calculated ana- 
lytically, and in terms of variables of — Xi±v are: 



Iq(v,w,x ,xi) = i ^2 + ^) (erf [of] - erf [a£ ]) + (°o ex Ph a o 2 ] ~ a i cxphaf 2 ]) , 

^(w.iu.aro.ai) =t\ (z + Z^j (erf[af] - erf [a$]) + exp[-a± 2 ] ^ T 3^ + ^ (l + 4') 

w 3 ( _|_ 



2^ 



± 2 i / w 
expl-ai 3-^3^] 



(7) 



(8) 
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I={v,w,x a , Xl ) = - 4 + 12 



(erf [of] - erf [4]) 



exp[ 



-a ±2 ]( T 8 



w _i_ 

12— of t: 



i + 4 2 ) + 



A u 



(3 + 2o ±2 )) 



1 . 4-2,/ W W _i_ U> / _i_2\ 10 4. / a 

Now, the collision frequency © can be easily calculated as 



±2 



I 



(9) 



(10) 



m=0 i=0 



We notice that this sum is made over all values of the 
index i, i.e. for all relative velocities x, but in practi- 
calapplications the situation is different. Integrals (JBJ) 
are homogenous functions of their arguments of the or- 
der zero, and we will now consider some of their general 
properties. Under conditions 



v > 6w, xq — 0, x\ = oo 
the integrals satisfy 

i+, it, 4 < io- 20 , 

Also, under conditions 

v > 6ui, Xo > u + 6w, X\ > Xq 
the integrals satisfy 

7 -, 7f, J 2 - <2-10- 16 . 



(11) 



(12) 
(13) 



(14) 



(15) 



Or in other words as v increases the integrals decrease 
very rapidly. 

When we work with a computer with double precision, 
machine accuracy is e m ~ 1CP 16 , so for a given ion ve- 
locity v > 6w, it is enough to integrate over the relative 
velocities up to v + 6w. Moreover, because of (fT2|) and 
(|13[) expansion ([3]) converges very rapidly. 

For ion velocities which are less than 6w similar 
consideration can be made. One can conclude, taking 
into account computer with a double precision, that the 
upper limit for integration in ([6]) should be at most 12w. 



A. Constant cross section 

For a constant cross section cr , from ©, © and ([7]) 
we obtain the collision frequency for the hard sphere scat- 



tering (which was previously found and used in |46|): 



v{v) = naow 
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: exp 



(f 



(16) 



From ([3]) we see that the collision frequency for a small 
ratio w/v (| 1 6 j) may be approximated by 



(17) 



u(v) = n<j v [1+2^2 



One can notice that, according to (J3j) , in the case of a con- 
stant cross section, only the correction of the first order 
exists. Relative difference between the last two formulas 
for the case v > 6w is less than machine accuracy e m 
and these formulas are represented by the same number 
in the computer register. 



B. Realistic cross sections 

Calculations of the collision frequency are also done 
for a real cross section (47| for Ar + ions in their parent 
gas. In Fig. [1] we show the collision frequency calculated 
by different methods. Exact solution and calculation by 
the equation [10] are essentially identical on the scale of 
the Figure. We can see that approximate formula from 
[43| should be improved further, as well as formula for 
the cold gas approximation. Also we can see that the 
analytic correction in ([3]), (which is included for v > 2w 
in Fig. |(T])) fits the exact formula very well. Relative dif- 
ference between the exact formula ^ and formula §5§ is 
less than 1.3 • 10 -3 for v > 2w, and less than 7.2 ■ 10 -6 
for v > 6w. Having in mind that relative difference be- 
tween the exact formula and the cold gas approximation 
@ even at v = 18w is greater than 10 -3 and increases 
with decreasing v we may conclude that corrected form 
of calculation is necessary. For large v/w limit and in 
order to save the computational time the analytic form 
([3"]) of the correction term is very useful and sufficiently 
accurate for all practical purposes. 
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FIG. 1: Collision frequency for Ar + ions in their parent gas 
for a cross sections from [471 ] for the parameter m = 2. Gas 
density was n = 10 22 m~ 3 and temperature T = 77 K. Solid 
line is calculated using exact formula (J2| or (110[) . dashed line 
using the formula from [4^], dash-dotted by using (|3} with 
inclusion of the correction terms after 2w, while the dotted 
curve is determined from the cold gas approximation Q. 



III. GAS VELOCITY DISTRIBUTION IN 
COLLISIONS WITH IONS 

When ion collides with gas, a special attention must be 
paid to statistical properties of gas particles which collide 
with that ion. Naively, one would think that statistical 
properties of gas particles which collide with an ion of 
a given velocity, are independent of that velocity and of 
the cross section. 

It has been customary in the literature 0, |48| to sam- 
ple gas velocities directly from the gas velocity distribu- 
tion function for their Monte Carlo simulations. This 
may lead to some inaccuracy especially for small ener- 
gies of particles. 

In that case they do not determine the distribution 
function of the collision probability (frequency) properly, 
as required in {T]). I n other words the time between two 
collisions is determined by the probability density for col- 
lisions which involves distribution function of the gas par- 
ticles but it should be modified to include the relative 
velocity dependence of the cross section. 

The error in sampling velocities from the gas velocity 
distribution function and applying them in the collision 
frequency may be small or may be compensated by other 
considerations but in that case one is never sure just how 
big the uncertainty is and whether it may become im- 
portant under some circumstances. We will see latter 
that such procedure will lead to results that are obvi- 
ously wrong (see Fig. 2]). 

While most authors use cold gas approximation and 
study mainly higher energies some authors have at- 
tempted to obtain a more accurate method of calculation 
of the distribution function of collision frequency in their 
calculations 0, S El • 

The probability density per unit time for collision of 



ion of velocity v with gas of velocity u is given by: 

P (v,u) = nfT(|v - u|) 'r u|F(u) , (18) 
v(v) 

where F(u) is the gas velocity distribution function and 
v{v) is the collision frequency, previously defined in (fTj). 
This is a conditional probability. We have ion of velocity 
v and we want velocity distribution of gas particles which 
collide with that ion. 

Randomly chosen gas velocity should be obtained from 
equation (fl8| . Method of gas velocity sampling from ([T8"|) 
depends on the shape of the cross section. One way to 
do it is to use the rejection method [49j which is easy 
for implementation, but often it is not efficient. Never- 
theless, it is very hard to develop an efficient and quite 
general method for sampling of a wide range of shapes of 
the cross sections so the rejection method is very often a 
method of choice. 

Two different implementations of the rejection method 
for (jl8p are obvious. Let q m ax be the maximum of 
<r(|v — u|) |v — u| for all allowed v and u (which are 
specified in the beginning of simulation) . This maximum 
can be chosen as 



V27 w 2 + 6vw + v 2 



(|v-u|). (19) 



Sampled velocity u' from the MVDF is accepted as the 
background velocity if 



cr(|v-u'|)|v-u'| > rq„ 



(20) 



is fulfilled (r is uniform random number from (0, 1]); oth- 
erwise we sample another velocity from the MVDF until 
the condition (120|) becomes fulfilled. This method is good 
when the cross section resembles a constant cross section 
and does not have sharp maxima. On the other hand, 
when the cross section has a sharp maximum for a very 
narrow interval of velocities (e. g. a(v) ~ previous 
method is very inefficient and we may use another esti- 
mation of q max simply as max VT .<T(v r )?; r , where v r goes 
over all allowed relative velocities (which are specified at 
the beginning of simulation). 



A. Constant cross section 

In the case of a constant cross section oo expression 
(dHJ) can be integrated analytically, and the probability 
density per unit time for collision of ion of velocity v with 
gas of velocity u is: 



p<(v,u) 



p>(v,u) 



na 



exp 



w \ UT 



(21) 



rtCTo 4w 



exp 



(22) 
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FIG. 2: Probability density ol the sampled gas velocity it: the 
case « = is shown by a dashed line, the case v — w is shown 
by a dotted line. MVDF is represented by a solid line. 



case when ion velocity is much higher than the thermal 
velocity of the background gas. In the case v = w we 
have the highest anisotropy and angles closer to tt are 
more probable. As ratio v/w increases or decreases an- 
gular distribution approaches the uniform case. Had we 
sampled the gas velocity from the gas velocity distribu- 
tion function, we would have obtained uniform azimuthal 
distribution, since the gas velocity distribution function 
is uniform in the space. 



IV. MIXTURES OF GASES 

If the background gas of density n is a mixture of k 
gases with relative densities on, 53;= i a i ~ 1j collision 
frequency Q becomes 




FIG. 3: Probability density of the azimuthal angle 9 in the 
case: v = - dashed line, v — w - dotted line, v = 5w or 
w = 5v - solid line, v = 2w or w = 2v - dash-dotted line. 



where p<(i>, u) means v < u and p>(i> , u) means v > u. 

In Fig. [5] we show the probability density of the sam- 
pled gas velocities for two different ion velocities in com- 
parison to the standard Maxwell's velocity distribution 
function. We see that, with an increasing ion velocity, 
we approach the MVDF. For ion velocities of the order 
of gas thermal velocity w and smaller, sampled gas ve- 
locities are considerably shifted towards higher velocities 
(in respect to the MVDF). 

Angular distribution may also be obtained from ex- 
pression (fl8|) . It is easy to see that the polar angle distri- 
bution is uniform, but the azimuthal angle distribution 
normalized by sinQ is nonuniform. Azimuthal angular 
dependence for a given ion velocity v and gas velocity V 



p(6) 



where cos 9 



3vV (v 2 + V 2 -2vVcos9) 



1/2 



+ v) 3 - \ v - 



(23) 



vV 

vV ■ 



In Fig. [3] we see that the uniform 
(isotropic) normalized azimuthal angle distributions cor- 
respond to the case of zero ion velocity and also to the 



v{v) = 2jaii/j(«), 



(24) 



Vi(v) =n \\ \ ^ (|v - u|) |v - u|Fi(u)du, (25) 



where Fj(u) is the gas velocity distribution (which has 
the same functional form for all gas components if the 
background gas is in equilibrium; the difference is only in 
the value of the thermal velocity due to different masses) 
and (7, is the ion-neutral cross section of the component i. 
Expression (|25[) can be easily calculated having in mind 
the previously described procedure for a pure background 
gas. 

Probability density per unit time for collisions of ions 
of velocity v with a gas particle of velocity u is given by 

p(y,u) S X: tt , W<(lv - U l ) Jr U|%) , (26) 



1=1 
v{v) 



Vi{v) 



(27) 



Sampling a uniform number in accordance to weighting 
factors Wi, we first determine the constituent of the mix- 
ture with which the ion of velocity v collides, and then 
we sample the velocity u of that constituent from the 
appropriate probability density 49]. It is a straightfor- 
ward procedure similar to that for a single component 
background gas. 



V. NUMERICAL RESULTS 

We have developed a new Monte Carlo simulation code 
for ion transport in background gas in external electric 
field. Since the Monte Carlo procedure is described else- 
where (see e.g. |50|) we will just give the part for colli- 
sion kinematics, because we have not found formulas in 
compact and simple form in literature valid for particle 
collision of arbitrary masses. 
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We consider two particles of masses m and M which 
collide. We assume that they have velocities v i V, re- 
spectively with respect to the laboratory coordinate sys- 
tem with e x , e y i e z . Assuming the momentum is con- 
served during collisions, the postcollision velocities are 



M 



where g' 

mM 



V 



9 e s 



m + M 
m 



m + M 



j(g-g'), 



(28) 
(29) 



V 



(v-V) 2 -fe j( , 



/' 



Mv+MV 

m+M 



,,.»,. , , , g = v — V, while E is in- 

elastic threshold. In the case of isotropic scattering e g > is 
isotropic unit vector which in simulation may be chosen 
using the algorithm proposed in [51] . 

In the case of anisotropic scattering let x (0 < X < 7r ) 
denote azimuthal and ip (0 < ip < 2tt) polar angle of 
the relative particle velocity after collision with respect 
to the velocity before collision. While polar angle is uni- 
form ip — 2nr, azimuthal is obtained as a solution of the 
equation 



Q°{ V -,X) sin xdx 
Jo a ( v 'X) sin xdx 



(30) 



where r is uniform random number from [0,1], and 
a(v,x) is the differential cross section. Then the post- 
collision velocities are 



V = V 



M ( g' . 
m + M ^g+-(hsmx-gcos X ) 

171 ( 5'/, • \ 

— T7 g+ — (hsmx- gcosx 

m + M \ g 



where h has components 



h x 
h y = 

K = 



9y + 9z sin ip, 
g x g y sinip + gg z cos ip 



g 2 y 



- n 



g x g z simp - gg y cos ip 



(31) 
(32) 

(33) 
(34) 

(35) 



9v+9 z 



Equivalent expressions for the case of purely elastic scat- 
tering (E = 0) can be found in [Ho| . 

The code was required to pass all the benchmark test 
calculations that are available for the zero gas tempera- 
ture. We implemented in the code the procedure for col- 
lision frequency calculation and for sampling the gas ve- 
locity as described in this paper. Number of benchmark 
results that could be used to test the non-zero gas tem- 
perature calculations that are available in the literature 
is very small. Thus we have tried to produce new values 
while at the same time repeating the available benchmark 
results. 

There are several sources for benchmark calculations 
for ions [I(J, but we found that only reference (53| 



does not use the cold gas approximation while providing 
a range of results that may be used to test the code. 

Numerical results of simulations in case of a constant 
cross section, for two ion-gas mass ratios are presented 
in Table [TJ There we compare our results to the Boltz- 
mann equation results of White et al. (53|. In case of 
a small ion-gas mass ratio and a constant cross section 
Monte Carlo simulation takes too long time because of 
very slow relaxation of ion transport properties (very in- 
efficient energy transfer from ion to the gas particle in 
elastic scattering because of the low mass ratio) , and our 
results are less accurate than the results for equal masses. 
Nevertheless, our results are in excellent agreement with 
the results from numerical solution of Boltzmann equa- 
tion [13 . 

In Fig. d] we show results from the simulation for dif- 
ferent electric fields. We applied different samplings pro- 
cedures for the background gas: from the correct for- 
mula (|18p and directly from the gas velocity distribu- 
tion function (which is assumed to be MVDF). One can 
see that thermal equilibrium (at zero field) cannot be 
achieved when the background velocities are not sam- 
pled in a proper way: mean ion energy is lower than the 
thermal energy |fcT and characteristic energy eDx/fJ. is 
lower than kT. 

In Fig. Owe show results for reduced mobilities for Ar + 
ions in their parent gas using the set of cross sections 
from [i^l • One can see that measured values at 293 K 
[54| are in excelent agreement with our calculated values, 
while measured values at 77 K [55[ are slightly lower than 
calculated. This may indicate the need to make some 
small adjustments to the cross section in the low energy 
limit. We should mention that measurements of Basurto 
et al. [56] agree very well with those of Helm and with 
our calculations. 

In Table [Til we tabulate the results for Ar + ions in 
their parent gas which may be used as benchmark results 
for verification of Monte Carlo codes or other techniques 
which take into account thermal effects of the background 
gas. We should mention that relative differences between 
results for two choices of parameters m — 2 and to = 2.3 
[47] ] are less that one percent, and only the results for 
to = 2 are presented here. 



VI. CONCLUSIONS 

We have developed an accurate and efficient algorithm 
for calculations of collision frequency in case when ther- 
mal motion of background gas cannot be neglected for 
the case of a Maxwellian velocity distribution of the back- 
ground gas particles. This is required for low energy sim- 
ulations of ion transport by using a Monte Carlo tech- 
nique or for implementation hybrid and in Particle in 
the Cell codes 0, [53] ■ Exact formulas are derived for 
very general cross section dependencies on relative ve- 
locity. Also, an analytic form of the correction terms 
in collision frequency was found, in the case of a large 
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TABLE I: Results for the ions with a hard sphere cross section model for m gas = 4amu, T = 293 K, E/n = 1 Td, <jq = 
6-10- 20 m 2 . 



e[eV] 

W [lO^s- 1 ] 
nD L [10 22 (ms) 
nD T [10 22 (ms) 
T T [10 2 K] 
T L [10 2 K] 



niton — nin 



Our results 



Results from [53] 



= 10" 



Our results 



Q Q, S 

Results from [531 



0.04272 
3.367 
0.885 
0.897 
3.074 
3.220 



0.04271 
3.368 
0.884° 
0.894° 
3.074 
3.220 



0.7332 
56.27 
155 
315 
56.70 
56.72 



0.73324 
56.187 
158.27 
313.20 
56.717 
56.727 



°Note an error of a factor 10 in the original data from the paper 



TABLE II: Benchmark results for Ar ions in Ar for m 







T = 77 K w 


= 179.0 ms" 1 




T = 293 K 


w = 349.2 ms -1 


E/n [Td] 


e [meV] 


v [ms — ] 


7iDL[10 22 (ms)^] 


e [meV] 


v [ids - ] 


uDl [10 22 (ms) -1 ] 


1 


9.97 


5.98 


00400 


37.89 


4.26 


0.0108 


5 


10.32 


29.86 


00410 


38.07 


21.22 


0.0108 


i n 








oo.u*± 


42.31 


U.U1UO 


16 


13.41 


92.07 


0.00482 


39.77 


67.21 


0.0110 


18 


14.23 


102.6 


0.00499 


40.26 


75.41 


0.0111 


20 


15.12 


113.0 


0.00516 


40.80 


83.48 


0.0111 


30 


20.25 


161.3 


0.00601 


44.13 


122.9 


0.0116 


35 


23.16 


183.5 


0.00644 


46.16 


141.9 


0.0119 


40 


26.24 


204.7 


0.00683 


48.38 


160.3 


0.0121 


45 


29.46 


224.8 


0.00724 


50.78 


178.3 
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ratio of ion - gas thermal velocities. It was emphasized 
that sampling of the gas velocity for a given ion veloc- 
ity (in Monte Carlo simulation) should be done by using 
appropriate probability of collisions. Explicit analytic 
formulas are calculated in the case of a constant cross 
section, and it was shown that the effective velocity dis- 
tribution function of the gas particles, which collide with 
an ion of a given velocity, is significantly different from a 
Maxwcllian. Results of our Monte Carlo simulation for 
ion transport are in excellent agreement with the numer- 
ical solution of Boltzmann equation and may be used as 
a benchmarks for ion transport properties. 

In recent literature there has been a greater tendency 
to represent ionic motion by a more complex model than 
the basic Langevin theory. Nanbu and coworkers have de- 



fined theory in terms of collision probabilities [58, 59, 60]. 
It may be worth the effort to include the presently pro- 
posed technique together with their approach, especially 
since that theory is normalized by using the mobilities 
which are usually available only in the thermal limit or 
close to it. 

It is worth noting that so far very accurate numerical 
techniques for solutions to the Boltzmann equation have 
been used to study ion transport [fill l62l . [fj3j . Applica- 
tion of Monte Carlo simulations for ions have not been 
as widespread for ions as it was for electrons. This to 
some degree is due to the problems in representing ther- 
mal collisions which are important over a much broader 
range of E/n. The present technique has already been 
applied in simulations of the transport of thermal ions 
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FIG. 4: Mean ion energy, drift velocity, longitudinal and transverse diffusion coefficient and transverse diffusion coefficient over 
mobility as a function of reduced electric field for a hard sphere model for rrii 0n — m gas = 4amu, T = 293 K, uq = 6 ■ 10 -20 m 2 
for different sampling procedures of velocity of background gas: from (|18[) - solid line; from MVDF - dashed line. 
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FIG. 5: Calculated reduced mobilities of Ar + ions in Ar at 
77 K (upper curve) and 293 K (lower curve). By dots are given 
experimental results from [55|. 



energies and values of E/n. Thus it seems better to apply 
Monte Carlo simulations to convert the transport data to 
cross sections as such codes may be directly compatible 
with plasma models or even more so directly applicable 
to build hybrid models that treat ion transport by MC 
simulations [Il,[33,[3l,[39|,[65, 66]. The present technique 
enables such application at a high level of accuracy with- 
out worrying about the effecs of thermal collisions on the 
determination of collision frequencies. It is also foresee- 
able that a similar implementation would be required for 
fast neutrals [65j, l66j . Just how big the effects of the inad- 
equate sampling of collision frequency on plasma model- 
ing are it is still not certain. Yet the procedure is simple 
enough that if once implemented it may serve as the ba- 
sis for hybrid or kinetic codes (MC section of PIC codes) 
that would leave little room to worry about the imple- 
mentation of proper technique and at the same time it 
would not burden the code. 



[36L f&3 | . At the same time the Boltzmann techniques 
mostly use interaction potentials to communicate the re- 
sults to the scientific community. While it is possible to 
convert those to the cross sections it is difficult to stan- 
dard modelers of plasmas to appreciate the applicability 
of certain sets of initial data over a wide range of mean 
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APPENDIX A 



Consider an integral of the form 



1 



f(x) exp 



[x — v) 2 



dx. 



(Al) 



Let f(x) be a function which is continuous and repeat- 
edly differentiable at x = v, and furthermore the integral 
in (|A1[) converges. We transform this integral by the 
substitution t — : 

w 

I = — I f(v + wt) exp \~t 2 ] dt. (A2) 

V^J-oo 

After expanding the function f(v + wt) into Taylor series 
f( v + wt ) = jr,^ffW(v), (A3) 
using the formula 

/ t 2n exp [-t 2 ] dx = r(n + -), (A4) 

J-oo * 



we obtain the final formula 



n=0 



(A5) 



In the limit of vanishing iu, the integral (|A1|) can be cal- 
culated directly using representation of delta function 



6(v) 



1 exp 
-= lim 

JTT w-tO 



(A6) 



When one wants to calculate an integral of the form 



V = 



/7TW 



f(x) exp 



(x — v) 2 



dx (A7) 



the situation is different. It is hard to obtain a general 
formula for integral but it is possible to calculate it 
in cases when the condition v ^> w holds. In that case 
main contribution to I' comes from points x which are 
in vicinity of v, so the lower limit in integral (|A7j) can be 
replaced by minus infinity, and I' is equal to the integral 
I up to the exponentially vanishing term: 



r = i + oi Exp 



(A8) 
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